---
title: "Information & Democracy, Replication File"
output: html_notebook
---

The script that follows replicates many of the analyses and graphics in Information and Democracy. It is distributed with the the following datasets:

Human Coding Datasets:

(1) all.coding.Rdata - a dataset including human coding of policy sentences

(2) all.pastvpresent.coding.Rdata - a dataset including past / present / future coding of policy sentences

Time Series Datasets:

(3) TSdata.Rdata (also distributed as TSdata.dta) - a time-series dataset of media, opinion and policy series by fiscal year

(4) TSdata.by.newsp.Rdata - a time-series dataset of individual newspaper signals

(5) TSdata.by.tv.Rdata - a time-series dataset of individual tv channel signals

Some of the analyses in Chapters 3 and 4 are based on large datasets of raw media data, and those analyses are not replicated below. This is because media databases are subject to licensing agreements and not publicly distributable in their original form. We also do not replicate analyses or distribute the Facebook data, which was provided to us by Dan Hiaeshutter-Rice. Even so, analyses that rely directly on either crowd-sourced human coding or on the time series derived from content analyses, polling and spending data are replicated here.

Note that all graphics and content analyses were produced in R, but most time series analysis was conducted in STATA. The script that follows uses the RStata package to produce results from STATA from within R, but the scripts for STATA can also be pasted into a STATA syntax file and run there using the TSdata.dta file.

Figures are not shown in this Markdown file, but are generated and will be saved in the folder alongside the data and this script.


```{r prep}

library(DataCombine)
library(quanteda)
library(stringr)
library(plyr)
library(tidyr)
library(stargazer)
library(RStata)
options("RStata.StataPath" = '/Applications/Stata/StataMP.app/Contents/MacOS/stata-mp')
options("RStata.StataVersion" = 15)

```


## Chapter 3


```{r Chapter 3 data}

load("all.coding.Rdata")

```


### Table 3.2
```{r Table 3.2}

#how many coded articles?
all <- c(length(C$uniqueid[C$domain=="def"]),
           length(C$uniqueid[C$domain=="wel"]),
           length(C$uniqueid[C$domain=="hea"]),
           length(C$uniqueid[C$domain=="edu"]),
           length(C$uniqueid[C$domain=="env"]))

#how many articles coded by Mturkers?
mturk <- c(length(C$uniqueid[C$domain=="def" & !is.na(C$relevance.m)]),
           length(C$uniqueid[C$domain=="wel" & !is.na(C$relevance.m)]),
           length(C$uniqueid[C$domain=="hea" & !is.na(C$relevance.m)]),
           length(C$uniqueid[C$domain=="edu" & !is.na(C$relevance.m)]),
           length(C$uniqueid[C$domain=="env" & !is.na(C$relevance.m)]))

#how many articles coded by students?
students <- c(length(C$uniqueid[C$domain=="def" & !is.na(C$relevance.s)]),
              length(C$uniqueid[C$domain=="wel" & !is.na(C$relevance.s)]),
              length(C$uniqueid[C$domain=="hea" & !is.na(C$relevance.s)]),
              length(C$uniqueid[C$domain=="edu" & !is.na(C$relevance.s)]),
              length(C$uniqueid[C$domain=="env" & !is.na(C$relevance.s)]))

tab <- as.data.frame(rbind(all,mturk, students))
colnames(tab) <- c("Defense","Welfare","Health","Education","Environment")
tab <- (t((tab)))
noquote(as.matrix(tab))

```

### Table 3.3
```{r Table 3.3 }

C$relevance.m.bin <- ifelse(C$relevance.m>.6,1,0)
C$dict.rel.bin <- ifelse(C$rel>0,"dict.yes","dict.no")
C$dict.dir.bin <- ifelse(C$up>0 | C$dn>0,1,0)

A <- table(C$relevance.m.bin[C$domain=="def"],C$dict.rel.bin[C$domain=="def"])
#class precision (number of correct predictions out of all predictions for that class)
def.precision.yes <- (A[2,2] / sum(A[,2])) 
#class recall (number of correct predictions out of number of documents humans coded in that class)
def.recall.yes <- (A[2,2] / sum(A[2,])) 
def.stat <- c(def.precision.yes,def.recall.yes)

A <- table(C$relevance.m.bin[C$domain=="wel"],C$dict.rel.bin[C$domain=="wel"])
wel.precision.yes <- (A[2,2] / sum(A[,2])) 
wel.recall.yes <- (A[2,2] / sum(A[2,])) 
wel.stat <- c(wel.precision.yes,wel.recall.yes)

A <- table(C$relevance.m.bin[C$domain=="hea"],C$dict.rel.bin[C$domain=="hea"])
hea.precision.yes <- (A[2,2] / sum(A[,2])) 
hea.recall.yes <- (A[2,2] / sum(A[2,])) 
hea.stat <- c(hea.precision.yes,hea.recall.yes)

A <- table(C$relevance.m.bin[C$domain=="edu"],C$dict.rel.bin[C$domain=="edu"])
edu.precision.yes <- (A[2,2] / sum(A[,2])) 
edu.recall.yes <- (A[2,2] / sum(A[2,])) 
edu.stat <- c(edu.precision.yes,edu.recall.yes)

A <- table(C$relevance.m.bin[C$domain=="env"],C$dict.rel.bin[C$domain=="env"])
env.precision.yes <- (A[2,2] / sum(A[,2])) 
env.recall.yes <- (A[2,2] / sum(A[2,])) 
env.stat <- c(env.precision.yes,env.recall.yes)

tab <- rbind(def.stat,wel.stat,hea.stat,edu.stat,env.stat)
colnames(tab) <- c("precision.yes","recall.yes")
rownames(tab) <- c("Defense","Welfare","Health","Education","Environment")
tab1 <- round(tab,2)

A <- table(C$relevance.m.bin[C$domain=="def" & C$dict.dir.bin==1],C$dict.rel.bin[C$domain=="def" & C$dict.dir.bin==1])
#class precision (number of correct predictions out of all predictions for that class)
def.precision.yes <- (A[2,2] / sum(A[,2])) 
#class recall (number of correct predictions out of number of documents humans coded in that class)
def.recall.yes <- (A[2,2] / sum(A[2,])) 
def.stat <- c(def.precision.yes,def.recall.yes)

A <- table(C$relevance.m.bin[C$domain=="wel" & C$dict.dir.bin==1],C$dict.rel.bin[C$domain=="wel" & C$dict.dir.bin==1])
wel.precision.yes <- (A[2,2] / sum(A[,2])) 
wel.recall.yes <- (A[2,2] / sum(A[2,])) 
wel.stat <- c(wel.precision.yes,wel.recall.yes)

A <- table(C$relevance.m.bin[C$domain=="hea" & C$dict.dir.bin==1],C$dict.rel.bin[C$domain=="hea" & C$dict.dir.bin==1])
hea.precision.yes <- (A[2,2] / sum(A[,2])) 
hea.recall.yes <- (A[2,2] / sum(A[2,])) 
hea.stat <- c(hea.precision.yes,hea.recall.yes)

A <- table(C$relevance.m.bin[C$domain=="edu" & C$dict.dir.bin==1],C$dict.rel.bin[C$domain=="edu" & C$dict.dir.bin==1])
edu.precision.yes <- (A[2,2] / sum(A[,2])) 
edu.recall.yes <- (A[2,2] / sum(A[2,])) 
edu.stat <- c(edu.precision.yes,edu.recall.yes)

A <- table(C$relevance.m.bin[C$domain=="env" & C$dict.dir.bin==1],C$dict.rel.bin[C$domain=="env" & C$dict.dir.bin==1])
env.precision.yes <- (A[2,2] / sum(A[,2])) 
env.recall.yes <- (A[2,2] / sum(A[2,])) 
env.stat <- c(env.precision.yes,env.recall.yes)

tab <- rbind(def.stat,wel.stat,hea.stat,edu.stat,env.stat)
colnames(tab) <- c("precision.yes","recall.yes")
rownames(tab) <- c("Defense","Welfare","Health","Education","Environment")
tab2 <- round(tab,2)

noquote(as.matrix(tab1))
noquote(as.matrix(tab2))

```


## Chapter 4


```{r Chapter 4 data}

load("all.coding.Rdata")
load("all.pastvpresent.coding.Rdata")

```


### Table 4.1
```{r Table 4.1}

C$relevance.s.bin <- ifelse(C$relevance.s>.6,1,0)
C$direction.m.cat[C$direction.m<(-.4)] <- (-1)
C$direction.m.cat[C$direction.m>.4] <- 1
C$direction.m.cat[C$direction.m>(-.401) & C$direction.m<.401] <- 0
C$direction.s.cat[C$direction.s<(-.4)] <- (-1)
C$direction.s.cat[C$direction.s>.4] <- 1
C$direction.s.cat[C$direction.s>(-.401) & C$direction.s<.401] <- 0

relevance.cor <- cor(C$relevance.s,C$relevance.m,use="complete.obs")
direction.cor <- cor(C$direction.s,C$direction.m,use="complete.obs")

A <- table(C$relevance.s.bin,C$relevance.m.bin)
agree <- A[1,1] + A[2,2] 
disagree <- A[2,1] + A[1,2] 
agree.p.all <- agree / (agree + disagree)
A <- table(C$relevance.s.bin[C$domain=="def"],C$relevance.m.bin[C$domain=="def"])
agree <- A[1,1] + A[2,2] 
disagree <- A[2,1] + A[1,2] 
agree.p.def <- agree / (agree + disagree)
A <- table(C$relevance.s.bin[C$domain=="wel"],C$relevance.m.bin[C$domain=="wel"])
agree <- A[1,1] + A[2,2] 
disagree <- A[2,1] + A[1,2] 
agree.p.wel <- agree / (agree + disagree)
A <- table(C$relevance.s.bin[C$domain=="hea"],C$relevance.m.bin[C$domain=="hea"])
agree <- A[1,1] + A[2,2] 
disagree <- A[2,1] + A[1,2] 
agree.p.hea <- agree / (agree + disagree)

A <- table(C$direction.s.cat,C$direction.m.cat)
agree <- A[1,1] + A[2,2] 
disagree <- A[2,1] + A[1,2] 
agree.p.all2 <- agree / (agree + disagree)
A <- table(C$direction.s.cat[C$domain=="def"],C$direction.m.cat[C$domain=="def"])
agree <- A[1,1] + A[2,2] 
disagree <- A[2,1] + A[1,2] 
agree.p.def2 <- agree / (agree + disagree)
A <- table(C$direction.s.cat[C$domain=="wel"],C$direction.m.cat[C$domain=="wel"])
agree <- A[1,1] + A[2,2] 
disagree <- A[2,1] + A[1,2] 
agree.p.wel2 <- agree / (agree + disagree)
A <- table(C$direction.s.cat[C$domain=="hea"],C$direction.m.cat[C$domain=="hea"])
agree <- A[1,1] + A[2,2] 
disagree <- A[2,1] + A[1,2] 
agree.p.hea2 <- agree / (agree + disagree)

df <- data.frame(matrix(vector(),5, 3,
                dimnames=list(c(), c("Domain", "Relevance","Direction"))),
                stringsAsFactors=F)
df$Domain <- c("All (cor)","All (% agreement)","Def (% agreement)","Wel (% agreement)","Hea (% agreement)")
df$Relevance <- round(c(relevance.cor,agree.p.all,agree.p.def,agree.p.wel,agree.p.hea),3)
df$Direction <- round(c(direction.cor,agree.p.all2,agree.p.def2,agree.p.wel2,agree.p.hea2),3)

noquote(as.matrix(df))

```


### Table 4.2
```{r Table 4.2}

A <- table(C$relevance.s.bin[C$domain=="def"],C$dict.rel.bin[C$domain=="def"])
#class precision (number of correct predictions out of all predictions for that class)
def.precision.yes <- (A[2,2] / sum(A[,2])) 
#class recall (number of correct predictions out of number of documents humans coded in that class)
def.recall.yes <- (A[2,2] / sum(A[2,])) 
def.stat <- c(def.precision.yes,def.recall.yes)

A <- table(C$relevance.s.bin[C$domain=="wel"],C$dict.rel.bin[C$domain=="wel"])
wel.precision.yes <- (A[2,2] / sum(A[,2])) 
wel.recall.yes <- (A[2,2] / sum(A[2,])) 
wel.stat <- c(wel.precision.yes,wel.recall.yes)

A <- table(C$relevance.s.bin[C$domain=="hea"],C$dict.rel.bin[C$domain=="hea"])
hea.precision.yes <- (A[2,2] / sum(A[,2])) 
hea.recall.yes <- (A[2,2] / sum(A[2,])) 
hea.stat <- c(hea.precision.yes,hea.recall.yes)

tab <- rbind(def.stat,wel.stat,hea.stat)
colnames(tab) <- c("precision.yes","recall.yes")
rownames(tab) <- c("Defense","Welfare","Health")
tab1 <- round(tab,2)

A <- table(C$relevance.m.bin[C$domain=="def" & !is.na(C$relevance.s.bin)],C$dict.rel.bin[C$domain=="def" & !is.na(C$relevance.s.bin)])
#class precision (number of correct predictions out of all predictions for that class)
def.precision.yes <- (A[2,2] / sum(A[,2])) 
#class recall (number of correct predictions out of number of documents humans coded in that class)
def.recall.yes <- (A[2,2] / sum(A[2,])) 
def.stat <- c(def.precision.yes,def.recall.yes)

A <- table(C$relevance.m.bin[C$domain=="wel" & !is.na(C$relevance.s.bin)],C$dict.rel.bin[C$domain=="wel" & !is.na(C$relevance.s.bin)])
wel.precision.yes <- (A[2,2] / sum(A[,2])) 
wel.recall.yes <- (A[2,2] / sum(A[2,])) 
wel.stat <- c(wel.precision.yes,wel.recall.yes)

A <- table(C$relevance.m.bin[C$domain=="hea" & !is.na(C$relevance.s.bin)],C$dict.rel.bin[C$domain=="hea" & !is.na(C$relevance.s.bin)])
hea.precision.yes <- (A[2,2] / sum(A[,2])) 
hea.recall.yes <- (A[2,2] / sum(A[2,])) 
hea.stat <- c(hea.precision.yes,hea.recall.yes)

tab <- rbind(def.stat,wel.stat,hea.stat)
colnames(tab) <- c("precision.yes","recall.yes")
rownames(tab) <- c("Defense","Welfare","Health")
tab2 <- round(tab,2)

noquote(as.matrix(tab1))
noquote(as.matrix(tab2))

```


### Figure 4.1
```{r Figure 4.1}

{
png("figure4.1.png",4,11, units = "in", res = 300)
par(mfrow=c(5,1))
par(mar = c(5.1, 4.1, 2, 2))
barplot(table(timing$timing[timing$domain=="def"])/length(timing$timing.c[timing$domain=="def"]),las=1,names.arg = c("Past","","<-","","","","","","->","","Future"),ylim=c(0,.2))
mtext("Defense",side=3,line=0,adj=0,cex=1.2)
mtext("Current",side=1,line=1,cex=.7)
par(mar = c(5.1, 4.1, 2, 2))
barplot(table(timing$timing[timing$domain=="wel"])/length(timing$timing.c[timing$domain=="wel"]),las=1,names.arg = c("Past","","<-","","","","","","->","","Future"),ylim=c(0,.2))
mtext("Welfare",side=3,line=0,adj=0,cex=1.2)
mtext("Current",side=1,line=1,cex=.7)
par(mar = c(5.1, 4.1, 2, 2))
barplot(table(timing$timing[timing$domain=="hea"])/length(timing$timing.c[timing$domain=="hea"]),las=1,names.arg = c("Past","","<-","","","","","","->","","Future"),ylim=c(0,.2))
mtext("Health",side=3,line=0,adj=0,cex=1.2)
mtext("Current",side=1,line=1,cex=.7)
par(mar = c(5.1, 4.1, 2, 2))
barplot(table(timing$timing[timing$domain=="edu"])/length(timing$timing.c[timing$domain=="edu"]),las=1,names.arg = c("Past","","<-","","","","","","->","","Future"),ylim=c(0,.2))
mtext("Education",side=3,line=0,adj=0,cex=1.2)
mtext("Current",side=1,line=1,cex=.7)
par(mar = c(5.1, 4.1, 2, 2))
barplot(table(timing$timing[timing$domain=="env"])/length(timing$timing.c[timing$domain=="env"]),las=1,names.arg = c("Past","","<-","","","","","","->","","Future"),ylim=c(0,.2))
mtext("Environment",side=3,line=0,adj=0,cex=1.2)
mtext("Current",side=1,line=1,cex=.7)
invisible(dev.off())
}

```


## Chapter 5


```{r Chapter 5 data}

load("TSdata.Rdata")
A$domain <- as.factor(A$domain)

A$Dsp2std <- A$sp2std - dplyr::lag(A$sp2std ,1)
A$sp2ch <- NA
A$sp2chstd <- NA
for (i in unique(A$domain)) {
  A$sp2ch[A$domain==i] <- A$sp2[A$domain==i] - dplyr::lag(A$sp2[A$domain==i])
  A$sp2chstd[A$domain==i] <- (A$sp2ch[A$domain==i] - mean(A$sp2ch[A$domain==i], na.rm=T)) / sd(A$sp2ch[A$domain==i], na.rm=T) 
}

load("TSdata.by.tv.Rdata")
load("TSdata.by.newsp.Rdata")

```


### Figure 5.1
```{r Figure 5.1}

{
pdf("figure5.1.pdf",7,12)
  
par(mfrow=c(5,1))

cor1 <- round(cor(A$sp2chstd[A$domain_text=="def" & !is.na(A$sig4n)],A$sig1n[A$domain_text=="def" & !is.na(A$sig4n)],use="complete.obs"),3)
cor2 <- round(cor(A$sp2chstd[A$domain_text=="def"],A$sig4n[A$domain_text=="def"],use="complete.obs"),3)
par(mar = c(5.1, 4.1, 2, 2.1))
plot(A$fy[A$fy>1979],A$sig1n[A$fy>1979],type="n",ann=F,axes=F,ylim=c(-3.5,3.5))
lines(A$fy[A$domain_text=="def"],A$sig1n[A$domain_text=="def"],col="gray",lwd=2,lty=1)
lines(A$fy[A$domain_text=="def"],A$sig4n[A$domain_text=="def"],col="gray",lwd=2,lty=2)
axis(1,col="gray",at=c(1980:2020))
axis(2,las=1,col="gray")
mtext("Media Signal / Spending Ch.",side=2,line=2,cex=.6)
mtext("Defense",side=3,line=0,adj=0,cex=1.1)
text(1980,3.5,"Correlation between spending change &",col="black",pos=4)
text(1982,2.8,paste0("Newspaper signal (solid): ",cor1),col="gray",pos=4)
text(1982,2.1,paste0("TV signal (dashed): ",cor2),col="gray",pos=4)
par(new=TRUE)
plot(A$fy[A$fy>1979],A$sp2chstd[A$fy>1979],type="n",ann=F,axes=F,ylim=c(-3.5,3.5))
lines(A$fy[A$domain_text=="def" & A$fy>1979],A$sp2chstd[A$domain_text=="def" & A$fy>1979],col="black",lwd=2)

cor1 <- round(cor(A$sp2chstd[A$domain_text=="wel" & !is.na(A$sig4n)],A$sig1n[A$domain_text=="wel" & !is.na(A$sig4n)],use="complete.obs"),3)
cor2 <- round(cor(A$sp2chstd[A$domain_text=="wel"],A$sig4n[A$domain_text=="wel"],use="complete.obs"),3)
par(mar = c(5.1, 4.1, 2, 2.1))
plot(A$fy[A$fy>1979],A$sig1n[A$fy>1979],type="n",ann=F,axes=F,ylim=c(-3.5,3.5))
lines(A$fy[A$domain_text=="wel"],A$sig1n[A$domain_text=="wel"],col="gray",lwd=2,lty=1)
lines(A$fy[A$domain_text=="wel"],A$sig4n[A$domain_text=="wel"],col="gray",lwd=2,lty=2)
axis(1,col="gray",at=c(1980:2020))
axis(2,las=1,col="gray")
mtext("Media Signal / Spending Ch.",side=2,line=2,cex=.6)
mtext("Welfare",side=3,line=0,adj=0,cex=1.1)
text(1980,3.5,"Correlation between spending change &",col="black",pos=4)
text(1982,2.8,paste0("Newspaper signal (solid): ",cor1),col="gray",pos=4)
text(1982,2.1,paste0("TV signal (dashed): ",cor2),col="gray",pos=4)
par(new=TRUE)
plot(A$fy[A$fy>1979],A$sp2chstd[A$fy>1979],type="n",ann=F,axes=F,ylim=c(-3.5,3.5))
lines(A$fy[A$domain_text=="wel" & A$fy>1979],A$sp2chstd[A$domain_text=="wel" & A$fy>1979],col="black",lwd=2)

cor1 <- round(cor(A$sp2chstd[A$domain_text=="hea" & !is.na(A$sig4n)],A$sig1n[A$domain_text=="hea" & !is.na(A$sig4n)],use="complete.obs"),3)
cor2 <- round(cor(A$sp2chstd[A$domain_text=="hea"],A$sig4n[A$domain_text=="hea"],use="complete.obs"),3)
par(mar = c(5.1, 4.1, 2, 2.1))
plot(A$fy[A$fy>1979],A$sig1n[A$fy>1979],type="n",ann=F,axes=F,ylim=c(-3.5,3.5))
lines(A$fy[A$domain_text=="hea"],A$sig1n[A$domain_text=="hea"],col="gray",lwd=2,lty=1)
lines(A$fy[A$domain_text=="hea"],A$sig4n[A$domain_text=="hea"],col="gray",lwd=2,lty=2)
axis(1,col="gray",at=c(1980:2020))
axis(2,las=1,col="gray")
mtext("Media Signal / Spending Ch.",side=2,line=2,cex=.6)
mtext("Health",side=3,line=0,adj=0,cex=1.1)
text(1980,3.5,"Correlation between spending change &",col="black",pos=4)
text(1982,2.8,paste0("Newspaper signal (solid): ",cor1),col="gray",pos=4)
text(1982,2.1,paste0("TV signal (dashed): ",cor2),col="gray",pos=4)
par(new=TRUE)
plot(A$fy[A$fy>1979],A$sp2chstd[A$fy>1979],type="n",ann=F,axes=F,ylim=c(-3.5,3.5))
lines(A$fy[A$domain_text=="hea" & A$fy>1979],A$sp2chstd[A$domain_text=="hea" & A$fy>1979],col="black",lwd=2)

cor1 <- round(cor(A$sp2chstd[A$domain_text=="edu" & !is.na(A$sig4n)],A$sig1n[A$domain_text=="edu" & !is.na(A$sig4n)],use="complete.obs"),3)
cor2 <- round(cor(A$sp2chstd[A$domain_text=="edu"],A$sig4n[A$domain_text=="edu"],use="complete.obs"),3)
par(mar = c(5.1, 4.1, 2, 2.1))
plot(A$fy[A$fy>1979],A$sig1n[A$fy>1979],type="n",ann=F,axes=F,ylim=c(-3.5,3.5))
lines(A$fy[A$domain_text=="edu"],A$sig1n[A$domain_text=="edu"],col="gray",lwd=2,lty=1)
lines(A$fy[A$domain_text=="edu"],A$sig4n[A$domain_text=="edu"],col="gray",lwd=2,lty=2)
axis(1,col="gray",at=c(1980:2020))
axis(2,las=1,col="gray")
mtext("Media Signal / Spending Ch.",side=2,line=2,cex=.6)
mtext("Education",side=3,line=0,adj=0,cex=1.1)
text(1980,3.5,"Correlation between spending change &",col="black",pos=4)
text(1982,2.8,paste0("Newspaper signal (solid): ",cor1),col="gray",pos=4)
text(1982,2.1,paste0("TV signal (dashed): ",cor2),col="gray",pos=4)
par(new=TRUE)
plot(A$fy[A$fy>1979],A$sp2chstd[A$fy>1979],type="n",ann=F,axes=F,ylim=c(-3.5,3.5))
lines(A$fy[A$domain_text=="edu" & A$fy>1979],A$sp2chstd[A$domain_text=="edu" & A$fy>1979],col="black",lwd=2)

cor1 <- round(cor(A$sp2chstd[A$domain_text=="env" & !is.na(A$sig4n)],A$sig1n[A$domain_text=="env" & !is.na(A$sig4n)],use="complete.obs"),3)
cor2 <- round(cor(A$sp2chstd[A$domain_text=="env"],A$sig4n[A$domain_text=="env"],use="complete.obs"),3)
par(mar = c(5.1, 4.1, 2, 2.1))
plot(A$fy[A$fy>1979],A$sig1n[A$fy>1979],type="n",ann=F,axes=F,ylim=c(-3.5,3.5))
lines(A$fy[A$domain_text=="env"],A$sig1n[A$domain_text=="env"],col="gray",lwd=2,lty=1)
lines(A$fy[A$domain_text=="env"],A$sig4n[A$domain_text=="env"],col="gray",lwd=2,lty=2)
axis(1,col="gray",at=c(1980:2020))
axis(2,las=1,col="gray")
mtext("Media Signal / Spending Ch.",side=2,line=2,cex=.6)
mtext("Environment",side=3,line=0,adj=0,cex=1.1)
text(1980,3.5,"Correlation between spending change &",col="black",pos=4)
text(1982,2.8,paste0("Newspaper signal (solid): ",cor1),col="gray",pos=4)
text(1982,2.1,paste0("TV signal (dashed): ",cor2),col="gray",pos=4)
par(new=TRUE)
plot(A$fy[A$fy>1979],A$sp2chstd[A$fy>1979],type="n",ann=F,axes=F,ylim=c(-3.5,3.5))
lines(A$fy[A$domain_text=="env" & A$fy>1979],A$sp2chstd[A$domain_text=="env" & A$fy>1979],col="black",lwd=2)
invisible(dev.off())
}

```


### Table 5.1
```{r Table 5.1}

stata_src <- "
    xtset domain fy
    quietly xtreg sig1n LD.sp2std, fe 
    estimates store B
    quietly xtreg sig1n D.sp2std, fe 
    estimates store C
    quietly xtreg sig1n FD.sp2std, fe 
    estimates store D
    esttab B C D, order(LD.sp2std D.sp2std FD.sp2std) cells(b(star fmt(3)) se(par fmt(3))) nogaps stat(r2 N)"
stata(stata_src,data.in=A)

```


### Table 5.2
```{r Table 5.2}

stata_src <- "
    xtset domain fy
    quietly xtregar sig1n L.D.sp2std D.sp2std F.D.sp2std, fe 
    estimates store E
    test L.D.sp2std = D.sp2std = F.D.sp2std
    quietly xtreg sig1n L.sig1n LD.sp2std D.sp2std FD.sp2std, fe 
    estimates store F    
    esttab E F, cells(b(star fmt(3)) se(par fmt(3))) nogaps stat(r2_w  N)"
stata(stata_src,data.in=A)

```


### Table 5.3
```{r Table 5.3}

stata_src <- "
    xtset domain fy
    quietly xtregar sig1n LD.sp2std D.sp2std FD.sp2std LD.sp1std D.sp1std FD.sp1std, fe 
    estimates store B
    quietly xtregar sig4n LD.sp2std D.sp2std FD.sp2std LD.sp1std D.sp1std FD.sp1std, fe 
    estimates store C
    esttab B C, cells(b(star fmt(3)) se(par fmt(3))) nogaps stat(r2_w N)"
stata(stata_src,data.in=A)

```


### Figure 5.2
```{r Figure 5.2}

df = data.frame(matrix(vector(), 0, 6,
                dimnames=list(c(), c("var", "coef", "stdev", "n", "rsq", "domain"))),
                stringsAsFactors=F)


for (w in 1:5) {
P <- A[A$domain==w,c("domain","fy","sig1n","sig2n","sig3n","sig4n","sp1std","sp2std")]
stata_src <- "
    tsset fy
    quietly reg sig1n LD.sp2std D.sp2std FD.sp2std
    regsave
    save results.dta, replace"
stata(stata_src,data.in=P, stata.echo=F)
temp <- rio::import("results.dta")
temp <- temp[,c("var","coef","stderr","N")]
temp$domain<-w
df <- rbind(df,temp)
}
np <- df

df = data.frame(matrix(vector(), 0, 6,
                dimnames=list(c(), c("var", "coef", "stdev", "n", "rsq", "domain"))),
                stringsAsFactors=F)
domains <- c("def","wel","hea","edu","env")  

for (w in 1:5) {
P <- A[A$domain==w,c("domain","fy","sig1n","sig2n","sig3n","sig4n","sp1std","sp2std")]
stata_src <- "
    tsset fy
    quietly reg sig4n LD.sp2std D.sp2std FD.sp2std
    regsave
    save results.dta, replace"
stata(stata_src,data.in=P, stata.echo=F)
temp <- rio::import("results.dta")
temp <- temp[,c("var","coef","stderr","N")]
temp$domain<-w
df <- rbind(df,temp)
}
tv <- df

{
pdf("figure5.2.pdf",6,15)
par(mfrow=c(5,2))

coefs <- np$coef[np$var!="_cons" & np$domain==1]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,3))
mtext("Defense",side=3,line=2,adj=0,cex=1.5)
mtext("Newspapers",side=3,line=(-1))
mtext("Change in Media Signal",side=2,line=2.5)

coefs <- tv$coef[tv$var!="_cons" & tv$domain==1]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,3))
mtext("TV",side=3,line=(-1))

coefs <- np$coef[np$var!="_cons" & np$domain==2]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,3))
mtext("Welfare",side=3,line=2,adj=0,cex=1.5)
mtext("Newspapers",side=3,line=(-1))
mtext("Change in Media Signal",side=2,line=2.5)

coefs <- tv$coef[tv$var!="_cons" & tv$domain==2]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,3))
mtext("TV",side=3,line=(-1))

coefs <- np$coef[np$var!="_cons" & np$domain==3]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,3))
mtext("Health",side=3,line=2,adj=0,cex=1.5)
mtext("Newspapers",side=3,line=(-1))
mtext("Change in Media Signal",side=2,line=2.5)

coefs <- tv$coef[tv$var!="_cons" & tv$domain==3]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,3))
mtext("TV",side=3,line=(-1))

coefs <- np$coef[np$var!="_cons" & np$domain==4]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,3))
mtext("Education",side=3,line=2,adj=0,cex=1.5)
mtext("Newspapers",side=3,line=(-1))
mtext("Change in Media Signal",side=2,line=2.5)

coefs <- tv$coef[tv$var!="_cons" & tv$domain==4]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,3))
mtext("TV",side=3,line=(-1))

coefs <- np$coef[np$var!="_cons" & np$domain==5]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,3))
mtext("Environment",side=3,line=2,adj=0,cex=1.5)
mtext("Newspapers",side=3,line=(-1))
mtext("Change in Media Signal",side=2,line=2.5)

coefs <- tv$coef[tv$var!="_cons" & tv$domain==5]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,3))
mtext("TV",side=3,line=(-1))

invisible(dev.off())
}

```


### Figure 5.3
```{r Figure 5.3}

C$domain[C$domain=="def"] <- "1"
C$domain[C$domain=="wel"] <- "2"
C$domain[C$domain=="hea"] <- "3"
C$domain[C$domain=="edu"] <- "4"
C$domain[C$domain=="env"] <- "5"
C$domain <- as.numeric(C$domain)
D$domain[D$domain=="def"] <- "1"
D$domain[D$domain=="wel"] <- "2"
D$domain[D$domain=="hea"] <- "3"
D$domain[D$domain=="edu"] <- "4"
D$domain[D$domain=="env"] <- "5"
D$domain <- as.numeric(D$domain)

E <- merge(A,C,by.x=c("fy","domain"),by.y=c("fy","domain"),all.x=T)
E <- merge(E,D,by.x=c("fy","domain"),by.y=c("fy","domain"),all.x=T)


E <- E[order(E$domain,E$fy),]
df = data.frame(matrix(vector(), 0, 7,
                dimnames=list(c(), c("var", "coef", "stdev", "n", "rsq", "domain","source"))),
                stringsAsFactors=F)
sources <- c("ajc","ark","azr","bgl","ctr","dvp","hch","lat","mst","nyt","ocr","phi","slp","stl","tbt","usa","wpo","abc","cbs","nbc","cnn","fox","msnbc")

for (w in 1:5) {
  for (i in sources) {
    d <- E[E$domain==w,c("fy",i,"sp2std")]
    colnames(d) <- c("fy","media","sp2std")
    d <- d[order(d$fy),]
    #d$media <- d$media/sd(d$media,na.rm=T)
    stata_src <- "
      tsset fy
      quietly prais media LD.sp2std D.sp2std FD.sp2std
      regsave
      save results.dta, replace"
    stata(stata_src,data.in=d,data.out=T, stata.echo=F)
    temp <- rio::import("results.dta")
    temp <- temp[,c("var","coef","stderr","N")]
    temp$domain<-w
    temp$source<-i
    df <- rbind(df,temp)
  }
}
rm(temp,d)

{
pdf("figure5.3.pdf",11,15)
par(mfrow=c(5,5))

coefs <- df$coef[df$var!="_cons" & df$domain==1 & df$source=="ajc"]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,500))
mtext("AJC",side=3,line=(-1))
mtext("Defense",side=3,line=2,adj=0,cex=1.5)
mtext("Coefficient",side=2,line=2.5)
coefs <- df$coef[df$var!="_cons" & df$domain==1 & df$source=="ctr"]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,500))
mtext("CTR",side=3,line=(-1))
coefs <- df$coef[df$var!="_cons" & df$domain==1 & df$source=="lat"]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,500))
mtext("LAT",side=3,line=(-1))
coefs <- df$coef[df$var!="_cons" & df$domain==1 & df$source=="nyt"]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,500))
mtext("NYT",side=3,line=(-1))
coefs <- df$coef[df$var!="_cons" & df$domain==1 & df$source=="wpo"]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,500))
mtext("WPO",side=3,line=(-1))

coefs <- df$coef[df$var!="_cons" & df$domain==2 & df$source=="ajc"]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,500))
mtext("AJC",side=3,line=(-1))
mtext("Welfare",side=3,line=2,adj=0,cex=1.5)
mtext("Coefficient",side=2,line=2.5)
coefs <- df$coef[df$var!="_cons" & df$domain==2 & df$source=="ctr"]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,500))
mtext("CTR",side=3,line=(-1))
coefs <- df$coef[df$var!="_cons" & df$domain==2 & df$source=="lat"]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,500))
mtext("LAT",side=3,line=(-1))
coefs <- df$coef[df$var!="_cons" & df$domain==2 & df$source=="nyt"]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,500))
mtext("NYT",side=3,line=(-1))
coefs <- df$coef[df$var!="_cons" & df$domain==2 & df$source=="wpo"]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,500))
mtext("WPO",side=3,line=(-1))

coefs <- df$coef[df$var!="_cons" & df$domain==3 & df$source=="ajc"]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,500))
mtext("AJC",side=3,line=(-1))
mtext("Health",side=3,line=2,adj=0,cex=1.5)
mtext("Coefficient",side=2,line=2.5)
coefs <- df$coef[df$var!="_cons" & df$domain==3 & df$source=="ctr"]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,500))
mtext("CTR",side=3,line=(-1))
coefs <- df$coef[df$var!="_cons" & df$domain==3 & df$source=="lat"]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,500))
mtext("LAT",side=3,line=(-1))
coefs <- df$coef[df$var!="_cons" & df$domain==3 & df$source=="nyt"]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,500))
mtext("NYT",side=3,line=(-1))
coefs <- df$coef[df$var!="_cons" & df$domain==3 & df$source=="wpo"]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,500))
mtext("WPO",side=3,line=(-1))

coefs <- df$coef[df$var!="_cons" & df$domain==4 & df$source=="ajc"]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,500))
mtext("AJC",side=3,line=(-1))
mtext("Education",side=3,line=2,adj=0,cex=1.5)
mtext("Coefficient",side=2,line=2.5)
coefs <- df$coef[df$var!="_cons" & df$domain==4 & df$source=="ctr"]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,500))
mtext("CTR",side=3,line=(-1))
coefs <- df$coef[df$var!="_cons" & df$domain==4 & df$source=="lat"]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,500))
mtext("LAT",side=3,line=(-1))
coefs <- df$coef[df$var!="_cons" & df$domain==4 & df$source=="nyt"]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,500))
mtext("NYT",side=3,line=(-1))
coefs <- df$coef[df$var!="_cons" & df$domain==4 & df$source=="wpo"]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,500))
mtext("WPO",side=3,line=(-1))

coefs <- df$coef[df$var!="_cons" & df$domain==5 & df$source=="ajc"]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,500))
mtext("AJC",side=3,line=(-1))
mtext("Environment",side=3,line=2,adj=0,cex=1.5)
mtext("Coefficient",side=2,line=2.5)
coefs <- df$coef[df$var!="_cons" & df$domain==5 & df$source=="ctr"]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,500))
mtext("CTR",side=3,line=(-1))
coefs <- df$coef[df$var!="_cons" & df$domain==5 & df$source=="lat"]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,500))
mtext("LAT",side=3,line=(-1))
coefs <- df$coef[df$var!="_cons" & df$domain==5 & df$source=="nyt"]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,500))
mtext("NYT",side=3,line=(-1))
coefs <- df$coef[df$var!="_cons" & df$domain==5 & df$source=="wpo"]
barplot(coefs,names.arg=c("t-1","t","t+1"),las=1,ylim=c(0,500))
mtext("WPO",side=3,line=(-1))

invisible(dev.off())
}

```


### Figure 5.4
```{r Figure 5.4}

df = data.frame(matrix(vector(), 0, 5,
                dimnames=list(c(), c("var", "coef", "stdev", "n", "source"))),
                stringsAsFactors=F)
sources <- c("ajc","ark","azr","bgl","ctr","dvp","hch","lat","mst","nyt","ocr","phi","slp","stl","tbt","usa","wpo","abc","cbs","nbc","cnn","fox","msnbc")

  for (i in sources) {
    d <- E[,c("fy","domain",i,"sp2std")]
    colnames(d) <- c("fy","domain","media","sp2std")
    d <- d[order(d$domain,d$fy),]
    for (w in 1:5){
      d$media[d$domain==w] <- d$media[d$domain==w]/sd(d$media[d$domain==w],na.rm=T)
    }
    stata_src <- "
      xtset domain fy
      gen LDsp2std = L.D.sp2std
      gen Dsp2std = D.sp2std
      gen FDsp2std = F.D.sp2std
      quietly xtregar media LDsp2std Dsp2std FDsp2std, fe
      lincom LDsp2std + Dsp2std + FDsp2std
      regsave
      save results.dta, replace"
    stata(stata_src,data.in=d,data.out=T, stata.echo=F)
    temp <- rio::import("results.dta")
    temp <- temp[,c("var","coef","stderr","N")]
    temp$source<-i
    df <- rbind(df,temp)
  }

rm(temp,d)

L <- read.csv("lincom.bysource.csv")

L <- L[order(L$coef),]
rownames(L) <- NULL
L$order <- as.numeric(rownames(L))

L$source <- stringr::str_to_upper(L$source)

{
pdf("figure5.4.pdf",6,8)
par(mar = c(5.1, 6.1, 2.1, 2.1))
plot(L$coef,L$order,type="n",ann=F,axes=F,xlim=c(-.5,2.5))
abline(v=0,lty=2,col="gray")
arrows(L$coef-(1.9*L$se),L$order,L$coef+(1.9*L$se),L$order,angle=90,length=.05,code=3,col="gray")
points(L$coef,L$order,pch=15,cex=1.5)
axis(1)
axis(2,las=1,at=c(1:23),labels=L$source)
mtext("Combined Impact of Spending (t-1, t, t+1),\nPooled Across All Policy Domains",side=1, line=3.5)
invisible(dev.off())
}


```


### Table 5.A1
```{r Table 5.A1}

stata_src <- "
    xtset domain fy
    gen pc4=d.sp2/l.sp2 *100
    cor sp2std pc4 if fy>1979
    quietly xtregar sig1n LD.sp2std D.sp2std FD.sp2std, fe
    estimates store A    
    quietly xtregar sig1n l.pc4 pc4 f.pc4, fe
    estimates store B
    quietly xtregar sig1n l.pc4 pc4 f.pc4 LD.sp2std D.sp2std FD.sp2std, fe
    estimates store C
    esttab A B C, cells(b(star fmt(3)) se(par fmt(3))) nogaps stat(r2_w  N) star(+ 0.10 * 0.05 ** 0.01 *** 0.001)"
stata(stata_src,data.in=A)


```


## Chapter 6


```{r Chapter 6 data}

load("TSdata.Rdata")

A$gini100b <- A$gini100
A$gini100b[A$domain==1] <- 0
A$gini100c <- A$gini100
A$gini100c[A$domain!=2] <- 0
A$rgdp <- A$gdp / A$deflator

for (i in 1:5) {
  A$cumsig1n[A$domain==i & A$fy>1979] <- cumsum(A$sig1n[A$domain==i & A$fy>1979])
  }
#lm(cumsig1n ~ sp2std, data=A[A$domain==1,])
#lm(cumsig1n ~ sp2std, data=A[A$domain==2,])
#lm(cumsig1n ~ sp2std, data=A[A$domain==3,])
#lm(cumsig1n ~ sp2std, data=A[A$domain==4,])
#lm(cumsig1n ~ sp2std, data=A[A$domain==5,])
A$Pcum <- NA
A$Pcum[A$domain==1] <- 3.390 * A$sp2std[A$domain==1]
A$Pcum[A$domain==2] <- 3.703 * A$sp2std[A$domain==2]
A$Pcum[A$domain==3] <- 15.057 * A$sp2std[A$domain==3]
A$Pcum[A$domain==4] <- 4.228 * A$sp2std[A$domain==4]
A$Pcum[A$domain==5] <- 17.139 * A$sp2std[A$domain==5]
A$Rcum <- A$cumsig1n - A$Pcum

```


### Table 6.1
```{r Table 6.1}

stata_src <- "
      xtset domain fy
      quietly xtreg pri rgdp gini100b sp2std, fe
      estimates store A
      quietly xtreg d.pri l.pri d.rgdp l.rgdp d.gini100b l.gini100b d.sp2std l.sp2std, fe
      estimates store B
      esttab A B, b(3) se(3) nogaps star(+ 0.10 * 0.05 ** 0.01 *** 0.001) r2 "
    stata(stata_src,data.in=A,data.out=T)

```


### Table 6.2
```{r Table 6.2}

stata_src <- "
      xtset domain fy
      quietly reg d.pri l.pri d.rgdp l.rgdp d.sp2std l.sp2std if domain==1
      estimates store A
      quietly reg d.pri l.pri d.rgdp l.rgdp d.sp2std l.sp2std d.gini100 l.gini100 if domain==2
      estimates store B
      quietly reg d.pri l.pri d.rgdp l.rgdp d.sp2std l.sp2std d.gini100 l.gini100 if domain==3
      estimates store C
      quietly reg d.pri l.pri d.rgdp l.rgdp d.sp2std l.sp2std d.gini100 l.gini100 if domain==4
      estimates store D
      quietly reg d.pri l.pri d.rgdp l.rgdp d.sp2std l.sp2std d.gini100 l.gini100 if domain==5
      estimates store E
      esttab A B C, b(3) se(3) nogaps star(+ 0.10 * 0.05 ** 0.01 *** 0.001) r2
      esttab D E, b(3) se(3) nogaps star(+ 0.10 * 0.05 ** 0.01 *** 0.001) r2"
    stata(stata_src,data.in=A,data.out=T)

```


### Table 6.3
```{r Table 6.3}

stata_src <- "
      xtset domain fy
      xtreg d.pri l.pri d.rgdp l.rgdp d.sp2std l.sp2std d.gini100c l.gini100c if fy>1980, fe
      estimates store A
      quietly xtreg d.pri l.pri d.rgdp l.rgdp d.sp2std l.sp2std d.gini100c l.gini100c sig1n l.cumsig1n, fe
      estimates store B
      quietly xtreg d.pri l.pri d.rgdp l.rgdp d.sp2std l.sp2std d.gini100c l.gini100c sig1n l.Pcum l.Rcum, fe
      estimates store C
      esttab A B C, b(3) se(3) nogaps star(+ 0.10 * 0.05 ** 0.01 *** 0.001) r2"
    stata(stata_src,data.in=A,data.out=T)

```


### Table 6.4
```{r Table 6.4}

stata_src <- "
      xtset domain fy
      quietly reg d.pri l.pri d.rgdp l.rgdp d.sp2std sig1n l.Pcum l.Rcum if domain==1 
      estimates store A
      quietly reg d.pri l.pri d.rgdp l.rgdp d.sp2std sig1n l.Pcum l.Rcum d.gini100 l.gini100 if domain==2
      estimates store B
      quietly reg d.pri l.pri d.rgdp l.rgdp d.sp2std sig1n l.Pcum l.Rcum if domain==3
      estimates store C
      quietly reg d.pri l.pri d.rgdp l.rgdp d.sp2std sig1n l.Pcum l.Rcum if domain==4
      estimates store D
      quietly reg d.pri l.pri d.rgdp l.rgdp d.sp2std sig1n l.Pcum l.Rcum if domain==5 
      estimates store E
      esttab A B C, b(3) se(3) nogaps star(+ 0.10 * 0.05 ** 0.01 *** 0.001) r2
      esttab D E, b(3) se(3) nogaps star(+ 0.10 * 0.05 ** 0.01 *** 0.001) r2"
stata(stata_src,data.in=A,data.out=T)

```


### Table 6.A1
```{r Table 6.A1}

stata_src <- "
      xtset domain fy
      gen pres = 0
      quietly replace pres =1 if fy==1980
      quietly replace pres =1 if fy>1992 & fy < 2001
      quietly replace pres =1 if fy>2008 & fy < 2017
      quietly replace pres = .5 if domain==1
      quietly xtreg d.pri l.pri d.rgdp l.rgdp d.sp2std l.sp2std d.gini100c l.gini100c sig1n l.cumsig1n d.pres l.pres, fe
      estimates store B
      quietly xtreg d.pri l.pri d.rgdp l.rgdp d.sp2std l.sp2std d.gini100c l.gini100c d.pres l.pres if e(sample), fe
      estimates store A
      esttab A B, b(3) se(3) nogaps star(+ 0.10 * 0.05 ** 0.01 *** 0.001) r2"
    stata(stata_src,data.in=A,data.out=T)

```


## Chapter 7


```{r Chapter 7 data}

load("TSdata.Rdata")

A$gini100b <- A$gini100
A$gini100b[A$domain==1] <- 0
A$gini100c <- A$gini100
A$gini100c[A$domain!=2] <- 0
A$rgdp <- A$gdp / A$deflator
A$circulation_th <-  (A$circulation *.001) -28.554

for (i in 1:5) {
  A$cumsig1n[A$domain==i & A$fy>1979] <- cumsum(A$sig1n[A$domain==i & A$fy>1979])
  }
#lm(cumsig1n ~ sp2std, data=A[A$domain==1,])
#lm(cumsig1n ~ sp2std, data=A[A$domain==2,])
#lm(cumsig1n ~ sp2std, data=A[A$domain==3,])
#lm(cumsig1n ~ sp2std, data=A[A$domain==4,])
#lm(cumsig1n ~ sp2std, data=A[A$domain==5,])
A$Pcum <- NA
A$Pcum[A$domain==1] <- 3.390 * A$sp2std[A$domain==1]
A$Pcum[A$domain==2] <- 3.703 * A$sp2std[A$domain==2]
A$Pcum[A$domain==3] <- 15.057 * A$sp2std[A$domain==3]
A$Pcum[A$domain==4] <- 4.228 * A$sp2std[A$domain==4]
A$Pcum[A$domain==5] <- 17.139 * A$sp2std[A$domain==5]
A$Rcum <- A$cumsig1n - A$Pcum

A$circcumsig1n <- A$circulation_th * A$cumsig1n

```


### Table 7.1
```{r Table 7.1}

stata_src <- "
      xtset domain fy
      quietly xtreg d.pri d.rgdp l.rgdp d.gini100c l.gini100c d.sp2std l.sp2std l.pri sig1n l.circulation_th  l.cumsig1n  l.circcumsig1n if domain<4, fe
      estimates store A
      esttab A, b(3) se(3) nogaps star(+ 0.10 * 0.05 ** 0.01 *** 0.001 ) r2"
stata(stata_src,data.in=A,data.out=T)

```


### Table 7.2
```{r Table 7.2}

stata_src <- "
      xtset domain fy
      gen dv=.
      replace dv = pr_np_lowi
      quietly xtreg d.dv d.rgdp l.rgdp d.gini100c l.gini100c d.sp2std l.sp2std sig1n l.cumsig1n L.dv if domain<4, fe
      estimates store A
      replace dv = pr_np_highi
      quietly xtreg d.dv d.rgdp l.rgdp d.gini100c l.gini100c d.sp2std l.sp2std sig1n l.cumsig1n L.dv if domain<4, fe
      estimates store B
      esttab A B, b(3) se(3) nogaps star(+ 0.10 * 0.05 ** 0.01 *** 0.001) r2"
stata(stata_src,data.in=A,data.out=T)

```


### Table 7.3
```{r Table 7.3}

stata_src <- "
      xtset domain fy
      gen dv = .
      replace dv = pr_ed_lowi
      quietly xtreg d.dv d.rgdp l.rgdp d.gini100c l.gini100c d.sp2std l.sp2std sig1n l.cumsig1n L.dv  if domain<4, fe
      estimates store A
      replace dv = pr_ed_highi
      quietly xtreg d.dv d.rgdp l.rgdp d.gini100c l.gini100c d.sp2std l.sp2std sig1n l.cumsig1n L.dv  if domain<4, fe
      estimates store B
      esttab A B, b(3) se(3) nogaps star(+ 0.10 * 0.05 ** 0.01 *** 0.001) r2"
stata(stata_src,data.in=A,data.out=T)

```


### Table 7.4
```{r Table 7.4}

stata_src <- "
      xtset domain fy
      gen dv = .
      replace dv = pr_pid_Di
      quietly xtreg d.dv d.rgdp l.rgdp d.gini100c l.gini100c d.sp2std l.sp2std sig1n l.cumsig1n L.dv  if domain<4, fe
      estimates store A
      replace dv = pr_pid_Ii
      quietly xtreg d.dv d.rgdp l.rgdp d.gini100c l.gini100c d.sp2std l.sp2std sig1n l.cumsig1n L.dv  if domain<4, fe
      estimates store B
      replace dv = pr_pid_Ri
      quietly xtreg d.dv d.rgdp l.rgdp d.gini100c l.gini100c d.sp2std l.sp2std sig1n l.cumsig1n L.dv  if domain<4, fe
      estimates store C
      esttab A B C, b(3) se(3) nogaps star(+ 0.10 * 0.05 ** 0.01 *** 0.001) r2"
stata(stata_src,data.in=A,data.out=T)

```


### Table 7.5-6
```{r Tables 7.5 and 7.6}

stata_src <- "
      xtset domain fy
      quietly regr sig1n l.sig1n d.sp2std i.domain if domain<4
      estimates store A
      quietly regr d.sp2std l2.sig1n l.d.sp2std i.domain  if domain<4
      estimates store B
      quietly regr pri l.pri l.sig1n i.domain  if domain<4
      estimates store C
      quietly regr sig1n l.sig1n l.pri i.domain  if domain<4
      estimates store D
      esttab A B, b(3) se(3) nogaps star(+ 0.10 * 0.05 ** 0.01 *** 0.001) r2
      esttab C D, b(3) se(3) nogaps star(+ 0.10 * 0.05 ** 0.01 *** 0.001) r2"
stata(stata_src,data.in=A,data.out=T)


stata_src <- "
      xtset domain fy
      quietly xtreg sig1n l.sig1n d.sp2std if domain<4, fe
      estimates store A
      quietly xtreg d.sp2std l2.sig1n l.d.sp2std  if domain<4, fe
      estimates store B
      quietly xtreg pri l.pri l.sig1n  if domain<4, fe
      estimates store C
      quietly xtreg sig1n l.sig1n l.pri  if domain<4, fe
      estimates store D
      esttab A B, b(3) se(3) nogaps star(+ 0.10 * 0.05 ** 0.01 *** 0.001) r2
      esttab C D, b(3) se(3) nogaps star(+ 0.10 * 0.05 ** 0.01 *** 0.001) r2"
stata(stata_src,data.in=A,data.out=T)

```

